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I. Abstract 

A recently developed numerical method based on the strip yield analysis approach was used to calculate the Crack 
Tip Opening Displacement (CTOD) and Crack Tip Opening Angle (CTOA) for a number of complex crack 
configurations of practical interest. This method is an adaptation of the dislocation-density based boundary element 
method. In this Boundary Element approach, crack-face opening displacements are obtained at any point on a crack 
using a series of definite integrals evaluated exactly in closed form for a variety of crack geometries, including 
infinite or finite extent, with arbitrarily applied loading conditions. 
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II. Introduction 

The Crack Tip Opening Displacement (CTOD) or, equivalently, the Crack Tip Opening Angle (CTOA), has become 
a reliable and convenient fracture criterion with which to assess crack extension and directionality in materials in 
which the effects of the plastic zone at the crack tip must be accounted for, e.g. in materials of considerable ductility 
and toughness. Many materials have been found 1 to exhibit a constant critical CTOD (or CTOA) from crack 
initiation through failure, making CTOD (or CTOA) a useful and reliable fracture criterion. 

Linear Elastic Fracture Mechanics (LEFM) was extended to materials in which significant plastic yielding occurred 
by Wells 2 , who postulated that CTOD was the variable which governs crack extension. Wells further proposed that 
the CTOD was proportional to the overall tensile strain even after general yield was reached. This proposition was 
later verified to be true, and became the basis of the widespread use of the CTOD in a variety of elastic -plastic 
assessments. 

The need for an explicit expression for the CTOD was filled when Dugdale 3 developed the strip yield model. 
Dugdale proposed that for a thin sheet loaded in tension, yielding would be limited to a narrow strip lying on the 
extension of the crack line, and that the effect of the cohesive yield loading would be to neutralize the stress 
singularity arising from the remote loading. Using Muskhelishvili 4 methods of complex-variable analysis, Dugdale’s 
method of strip yield analysis obtained an explicit expression for the CTOD, while circumventing the need to 
perform full-blown elastic-plastic analysis. It thus offers an efficient approach to modeling elastic-plastic behavior 
through the superimposition of two elastic solutions. 

Dugdale’s original development of the strip yield model involved a thin plate of infinite extent under plane stress 
conditions. While closed-form solutions have been obtained for other infinite geometries, these have necessarily 
been limited to a few particular cases; and efforts to extend the method to a general-purpose technique for arbitrary 
finite geometries have been hampered because the needed elastic solutions, in general, can be obtained only through 
painstaking numerical computations that have to be tailor-made for the problem at hand. 

While methods like the Weight function, Green’s function and Collocation methods 5 " 9 have been used to develop 
solutions for finite geometries, these have been techniques developed for particular configurations that also do not 
lend themselves to general-purpose use. Weight function methods, for instance, while capable of extension to 
arbitrary loading situations, require the calculation of two or more reference solutions, and impose thereby a heavy 
computational burden, often unacceptable to the practicing fracture analyst. 
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Finite-element analysis, while certainly general, also exacts a severe computational toll from the user: a full strip 
yield analysis requires determination of stress intensity for the actual (arbitrary) loading and the cohesive yield stress 
loading near the crack tip, together with the determination of the crack- face opening displacement at the junction of 
the load-free and yield-loaded parts. In addition, one must choose from several locations behind the crack tip which 
have been proposed for measuring CTOD: the first node 10 ' 12 , the second node 1213 , or the 45° intercept 14 . Finally, 
finite element modeling, especially when done for a range of crack lengths, imposes a severe re-meshing burden 
because of the multitude of meshes required, and is thus not readily adaptable to general-purpose strip yield analysis. 

To overcome the aforementioned shortcomings, a computational method has been developed 15 to obtain crack-face 
opening displacements at any point on a crack face, and thus has the capability to use the strip-yield model for 
general two-dimensional crack problems. This technique adapts the dislocation-density based boundary-element 
method (BEM) developed by Chang and Mear 16 in which the crack line integrals are treated in a special way to 
eliminate the difficulties associated with the application of the standard BEM to fracture problems. The resulting 
integral equations contain integrals associated with the fracture in terms of a distribution of point dislocations along 
the crack line. With this adaptation, presently implemented in the NASGRO® fracture mechanics analysis 
software 17 , the CTOD can be easily and accurately computed for arbitrary plane geometries and loadings, and the 
modeling requirements are free of the onerous modeling burden of other methods. 

This paper makes use of this novel approach to calculate the CTOD for several complex crack configurations of 
practical interest. 


III. Theory 

A direct application of conventional BEM to fracture problems leads to a mathematically degenerate formulation 
due to the geometric proximity of the crack surfaces, and information about the tractions on the crack is lost 18 . This 
can be circumvented by developing an integral equation for the crack face tractions, for example, by using the 
displacement integral equation for interior points ( i.e . Somigliana’s identity) and applying the strain-displacement 
relations, Flooke’s Law, and a proper limiting process as the interior point approaches the crack surface. Flowever, 
the resulting relation contains a hypersingular kernel which requires a special interpretation and is challenging to 
evaluate numerically. In addition, this formulation requires that the gradients of the relative crack face 
displacements be continuous, a difficult and computationally intensive requirement to implement 10 ). 

In an alternative approach proposed by Chang and Mear 16 , the hypersingularity is avoided by eliminating the 
Cauchy-type kernel in the displacement integral equation before forming the displacement gradients. This is 
accomplished in by an appropriate choice of stress function for the stress fundamental solution due to a point load, 
and an integration by parts to obtain a ‘■‘modified” displacement integral equation. Also, the integrals over the upper 
and lower crack surfaces can be written as integrals over a single crack surface, e.g. the upper surface. This results 
in the crack face displacement terms being replaced by the relative crack face displacements AD, where AD is the 
difference between the displacements of the upper and lower surfaces, i.e. the crack-opening displacement. A 
traction relation for points on the crack faces is then easily derived from this modified displacement equation in the 
usual manner. 

It can readily be shown that the boundary integral equations upon which this BEM is based consist of the 
conventional ordinary boundary terms (i.e. conventional displacement and traction fundamental solutions) and two 
additional terms which can be identified as distributions of point loads and of point dislocations along the crack 
lines. 

In this BEM the gradients of the relative crack face displacements are described by the dislocation density function 
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where the complex-valued relative displacement function is defined as 
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In (2), AD is the difference in displacements between the upper and lower crack surfaces, u x and u y are the x- and y- 
components, respectively, of displacement, K equals 3-4v for plane strain and (3-v)/(l+v) for plane stress, u and v 
are the shear modulus and Poisson’s ratio, respectively, and s is the distance along the crack from the tip to the point 
of interest. 

Integrating (2) yields nominally 
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where evaluation of the integral requires the summing up of the contribution from each element between the crack 
tip and the point of interest. 

In Ref. 16, the geometric interpolation along the crack line is piecewise linear; the crack line geometry over the 
(j+l)ih element is approximated by 
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The dislocation density function on each element is approximated by functions that contain the requisite singularity 
multiplied by an interpolating part that depends on the natural element coordinate t, 
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where it is assumed that the number of elements is even. In (5) and (6), Aj are nodal quantities, which are obtained 
from the vector of unknowns; a,j is the half-length of the yth element, and pj = / ?, /a, in which [ 3j is the element-wise 
arc length from the center of the /th element to the nearest crack tip. To ensure continuity of the dislocation density, 
the total arc length of the “left” portion of the crack must be equal to the total arc length of the “right” portion: 

Pj + OL j = p j+l + a j+l for j = 'AN (7) 

Substituting (4), (5) and (6) in (3) allows evaluation of AD as an element-wise sum of definite integrals. Most of the 
integrals thus obtained involve integrands that are non-singular and well behaved, and are generically of the form 
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These are integrated exactly in closed form using elementary methods. The significant point is that even those 
integrals which have the crack tips at the end-points of integration (and consequently integrands that are singular) 
have been integrated exactly in closed form, yielding results that are independent of the integration technique 
employed. 

Furthermore, the crack tip stress intensity factor (SIF) is directly related to the limiting value of the dislocation 
density function at the crack tip: 
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where 1 is the angle of inclination of the unit tangent vector at the crack tip. For the linear crack discretization 
described above, (10) simplifies to 
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for the left crack tip, where Aj is the nodal value of A(t) at the left crack tip. A similar expression can be written for 
the right crack tip. 

This approach can be applied to any crack geometry, of either infinite or finite extent, with arbitrary applied loading 
conditions, including mixed mode crack problems. 


IV. Procedure 


Dugdale proposed that all plastic deformation is concentrated in a strip in front of the crack and that this problem 
can be solved by superposition of two elastic solutions. The problem of a body with a crack of physical length 2a 
and plastic zone size p at each crack tip is equivalent to the same body with a crack of physical length 2(a+ p) where 
the portion of the crack flanks corresponding the plastic zones (of length p) are subject to a closure stress equal to 
the yield stress, as shown in Figure l 20 . 
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Figure 1. The Dugdale approach 


The strip-yield model then is a superposition of the two elastic solutions in plates C and D wherein the cohesive 
yield loading on the crack flanks neutralizes the stress singularity arising from the remote loading. 

The BEM analysis consists of constructing a model as shown in plate B with p chosen so that the SIF at the notional 
crack tip ( i.e . at a+ p) is 0. In practice, for a given yield stress, this can be achieved by either setting the plastic zone 
size and iterating on the remote stress, or setting the remote stress and iterating on the plastic zone size, until the SIF 
is zero. The fonner has the advantage of having to create a mesh for the crack and plastic zone just once, while the 
latter has the advantage of determining the CTOD for specific values of the remote stress chosen a priori. The 
CTOD value is then given by the crack face displacement at the intersection of the loaded and stress-free parts of the 
crack. 






The BEM formulation used herein uses quadratic boundary elements and linear crack elements; an additional 
advantage of this formulation is that it generally requires far fewer elements for highly accurate stress intensity 
factor, stress, and displacement analysis than other BEM formulations. For example, for the crack configurations 
discussed in this paper, accurate results with errors of 3% or less can be obtained typically by using fewer than 20 
elements each on the external and internal boundaries and crack line. The crack face loading discontinuity requires 
a finer mesh immediately to either side of the discontinuity than elsewhere on the crack, and it should be noted that 
in the absence of such a discontinuity even fewer crack elements would yield accurate results. 


V. Results 

In this section, several examples are presented to demonstrate the accuracy and versatility of this technique to 
calculate CTOD and CTOA. The examples include cases for finite and infinite domains, for center and edge cracks, 
and for complex configurations containing one or more holes. The verification cases are compared to well-known 
closed-form solutions as well as for experimental results for several common aerospace materials. 

In each example, the crack was meshed with linear elements in several segments of differing element densities in 
order to adequately capture the discontinuity of the crack face loading. External and hole boundaries were meshed 
with quadratic elements, with comparatively sparse density required for accuracy, as discussed in the previous 
section. 

A. Plastic zone size and CTOD for a crack in an infinite sheet 

The first verification case is a simulation of the original Dugdale strip-yield model, a crack in an infinite sheet 
subject to remote uniform tension with cohesive loading equal to the yield stress on the portions of the crack faces of 
length p near each tip as shown in Figure 2. 



Figure 2. The Dugdale model - a crack in an infinite sheet with cohesive crack face 

loading near the crack tips 


The physical crack is of length 2a while the notional crack to be modeled is of length 2a+2p. Results for a wide 
range of applied stress ratio are shown in Figure 3 for the CTOD (normalized by the yield stress a y , physical crack 
size a, and Y oung’s modulus E where E ’ is E for plane stress and E/(l-v 2 ) for plane strain) and for the relative 
crack size c/a, where c= a+p. In Dugdale ’s closed- form solution, the CTOD is given as 
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where the crack size ratio c/a is given by 
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Agreement between the BEM calculations and the Dugdale closed form solution is within 1%. 
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Figure 3. Comparison of the BEM to the Dugdale solution for relative crack size and normalized CTOD 


B. Plastic zone size for a middle crack in a finite-width sheet 

Figure 4 shows good agreement between NASBEM calculations and experimentally measured plastic zone sizes 21 
for 0.020” thick AM350CRT steel sheet. As a comparison, the Dugdale and Irwin plastic zone sizes are also shown. 
For this combination of width and crack size the finite-width correction factor differs from the infinite sheet 
solution’s value by less than 2%, so the small difference between the BEM calculations and the Dugdale infinite 
sheet solution is to be expected. Also, as is expected, the Dugdale and Irwin solutions bound the experimental 
results. 



Figure 4. Comparison of the BEM to experimentally measured plastic zone sizes 


C. Plastic zone size and CTOD for periodic cracks in an infinite sheet 

One of the long-term goals of the current research is to demonstrate the viability of using the Dugdale technique as 
implemented in NASGRO® to calculate CTOD and CTOA for use in characterizing the fracture behavior and 
linkup of multiple cracks in very wide panels such as fuselage panels. To this end, Tada’s solution for periodic 
crack in an infinite sheet 22 was chosen as a validation case. 



Figure 5. Configuration of the case of periodic cracks in an infinite sheet 


This case presents some unique issues of practical concern to the analyst, namely, how to model an infinite number 
of cracks? Guidance was taken from literature on stress intensity factor analysis of stiffened fuselage panels 23 in 
which it was found that a periodic set of stiffeners in an infinite sheet could be reasonably approximated by seven 
stiffeners. In this vein it is postulated that a periodic set of cracks in an infinite sheet can be adequately represented 
by 7 cracks. As a practical matter, however, due to program size limitations the outer two cracks on each side were 
combined to form a single crack of double size, so that this case was in fact modeled as 5 collinear cracks. As can 





be seen in Figures 6 and 7, the computational results match the closed- form results very well for a wide range of 
crack size to crack spacing ratios, both for computed plastic zone sizes as well as for CTOD. 
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Figure 6. Comparison of the BEM to Tada’s solution for periodic cracks in an infinite sheet 

for crack size to crack spacing ratios 


Comparison of BE to T ada - periodic cracks in infinite sheet 

CTOD/(o y *2b/E') 



c/b [actual crack + plastic zone / spacing] 


Figure 7. Comparison of the BEM to Tada’s solution for periodic cracks in an infinite sheet 

for normalized CTOD 


D. CTOA at onset of stable tearing for a middle crack in a finite-width sheet 

The second long-term goal of the current research is to use the Dugdale technique to predict fracture when there is 
large stable tearing by using the critical CTOA. CTOA is often seen as a better fracture criterion due to its local 
nature than global measures such is the /-integral, which can be specimen and crack size dependent. 

Figure 8 shows a typical distribution of CTOA values with respect to stable crack extension for A1 2024-T3 M(T) 
and C(T) specimens 24 . After an initial transient region (attributable to constraint and tunneling issues), the critical 
CTOA is constant over a large range of crack extension. It is this CTOA that is intended to be predicted in future 
research, and while one can calculate CTOA from CTOD, the Boundary Element code in NASGRO® in its present 
state is not yet able to simulate stable tearing. 



Figure 8. Typical distribution of experimentally measured CTOA as a function of stable crack extension 


In the preceding comparisons CTOD was taken as the crack opening displacement at the tip of the physical crack, 
i.e. at the intersection of the traction-free and the yield-loaded portions of the crack. As a practical matter, however, 
CTOD and CTOA are typically measured at some distance behind the physical crack tip, usually 0.02-0.06” (0.5- 
1.5mm). CTOA is calculated from CTOD by the well-known formula 


CTOA = 2tan 1 
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where d is the distance behind the crack tip at which CTOD is measured. 


Figure 9 shows that the CTOA calculated from the CTOD computed in the BEM at onset of stable tearing is 
reasonably constant for a range of initial crack sizes, and the values are similar to the data presented in Figure 8. 



Figure 9. BEM calculation of CTOA at onset of stable growth for 0.060” M(T) specimens 





E. CTOA at onset of stable tearing for a three-hole tension specimen 


The final case studied is that of a three-hole tension specimen. This specimen, shown in Figure 10, was designed to 
simulate cracks in stiffened panels, and Figure 1 1 shows the complexity of the stress intensity factor as a function of 
crack length. 



Figures 10 and 11. Configuration of the three-hole tension specimen and 
variation of the stress intensity factor with crack size 24 


Figure 12 shows the CTOA calculated by the BEM at onset of stable tearing for a number of initial crack sizes, and 
comparison to test results 24 in Figure 13 shows that the BEM values are within range. Also shown is the variation of 
failure load with initial crack size. 
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Figure 12. BEM calculation of CTOA at onset of stable growth 
and corresponding failure stresses for three-hole tension specimen 



Figure 13. Experimentally measured CTOA as a function of stable crack extension 
for a three-hole tension specimen" 4 


VI. Summary 

The feasibility of a new numerical method for calculating the Crack Tip Opening Displacement (CTOD) and Crack 
Tip Opening Angle (CTOA) was investigated. This new technique adapted the Boundary Element / Dislocation 
Density Method to obtain crack-face opening displacements at any point on the crack. Most of the methods in the 


literature used to calculate CTOD are complicated and resource intensive, and many of them apply only to specific 
geometries and loading conditions. The new approach developed in this paper was successfully applied to a number 
of different crack configurations including finite and infinite geometries containing complex geometrical features 
and arbitrary applied loading conditions. Results were obtained with reasonable accuracy, requiring a minimum of 
computational time and resources. 
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